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CHAPTER    1 
INTRODUCTION 

In  this  era  of  the  digital  computer,  scientists  and 
technologists  are  finding  new  ways  to  use  the  computer  as  a 
problem-solving  tool  in  their  areas.  This  is  especially 
true  of  the  engineering  sciences.  Engineering  problems 
involve  enormous,  somtimes  tedious,  calculations; 
therefore,  computers  are  exceptionally  well-suited  for 
applications     in    this    field. 

The  engineering  design  process  is  often  one  of  trial 
and  error;  however,  computer  simulation  has  done  much  to 
reduce  the  human  effort  involved  and  to  improve  efficiency. 
The  simulation  of  a  design  is  carried  out  with  computer 
software  --  the  so  called  simulator.  Different  simulators 
are  available  with  their  own  objectives.  In  general,  a 
simulator  is  a  machine  or  program  to  which  one  feeds  the 
description  of  the  system  and  obtains  the  response 
characteristics  of  the  system.  The  specific  simulator 
covered  in  this  research  is  one  for  simulation  of  robotic 
dynamics,  both  forward  and  inverse.  The  input  to  the 
simulator  will  be  the  robot  parameters  such  as  the  masses, 
moments  of  inertia,  and  dimensions  of  each  link,  joint 
types     and     locations,      etc.        For     the     forward     dynamic 


simulation,  it  also  requires  joint  torques  or  forces, 
initial  positions  and  velocities  of  each  link  to  find  the 
link  positions,  velocities,  and  accelerations  at  various 
instants  of  time.  In  the  inverse  dynamic  simulation,  link 
positions,  velocities,  and  accelerations  have  to  be 
specified;  as  the  simulation  proceeds,  it  will  find  the 
required  joint  torques  or  forces  at  each  time  interval. 

Due  to  difference  among  manipulator  designs,  simulators 
are  often  written  for  each  specific  design.  When  a 
designer  wants  to  change  the  basic  layout  of  his  design, 
the  simulator  also  needs  to  be  changed.  Even  minor 
modifications  in  the  robot  system  could  result  in 
significant  modification  of  the  simulator  program.  This  is 
felt  most  strongly  when  the  simulator  is  being  used  in  the 
design  of  a  robot,  since  it  is  likely  that  the  robot  will 
undergo  considerable  modification  during  the  design 
process.  This  lack  of  versatility  is  a  serious  handicap 
that  provides  the  motivation  for  the  development  of  more 
flexible  simulators.  A  good  way  to  do  this  job  is  to  write 
an  automatic  simulator  generator,  as  described  in  this 
thesis,  which  can  handle  a  wide  variety  of  robot 
manipulator  systems.  This  is  best  done  through  symbolic 
computation  using  an  algebraic  manipulation  language,  such 
as  MACSYMA  or  REDUCE-3. 

There  are  several  methods  available  for  setting  the 
dynamic  equations,   for  example.   Lagrange  [1,2,3],   virtual 


work  [4],  and  Newton-Euler  [5.6.7].  For  the  purpose  of 
control,  it  is  desirable  to  obtain  the  solution  in  real- 
time. To  solve  the  dynamic  equations  efficiently,  several 
ideas  have  been  presented  in  the  literature.  For  example, 
the  recursive  Lagrange  formulation  [8]  was  used  to  solve 
the  inverse  dynamic  problem,  in  which  the  generalized 
forces  were  expressed  in  recurrence  relations  involving  the 
generalized  coordinates,  velocities,  and  accelerations.  By 
the  recurrence  of  these  quantities,  it  is  possible  to 
reduce  the  number  of  multiplications  and  additions  involved 
from  quartic  to  linear  dependence  on  the  number  of  links  of 
the  system.  There  is  also  a  recursive  Newton-Euler 
formulation  [5,6,9],  which  uses  a  similar  idea.  Another 
approach  is  the  tabular  solution  look-up  [10].  which 
reduces  computation  but  increases  the  input/output  (I/O) 
requirements.  The  inverse  dynamics  problem  for  specific 
motions  of  a  specific  robot  is  solved  once  and  all  the 
solutions  are  stored,  then  one  can  find  the  solution  of 
generalized  forces  at  a  certain  instant  by  interpolation. 
For  specific  robots,  it  is  also  possible  to  simplify 
problems  by  dropping  Coriolis  and  centrifugal  terms  [1,11]. 
This  simplification  can  be  made  only  after  verifying  that 
these  terms  have  insignificant  contributions  and  that 
certain  presumptions  apply. 


Among  these  methods,  the  Lagrangian  formulation  was 
chosen  as  being  the  most  suitable  for  computer 
implementation.  One  coordinate  is  assigned  to  each  link, 
thus  we  are  dealing  with  the  minimum  set  of  generalized 
coordinates.  Trans lational  joints  and  revolute  joints  are 
modeled  in  a  unified  way;  this  streamlines  the  program 
implementation.  The  same  Lagrange  equation  formulation  is 
able  to  represent  both  forward  and  inverse  dynamic 
problems;  some  common  terms  have  the  same  expressions  in 
both  cases  which  in  turn  further  simplifies  the 
implementation.  The  forces  resulting  from  springs  and 
dampers,  frictional  forces,  and  other  externally  applied 
forces  are  taken  into  account  in  the  generalized  force 
term . 

As  stated  previously,  the  forward  and  inverse  dynamic 
problems  are  to  be  solved.  For  the  inverse  one,  the 
situation  becomes  straight  forward  since  all  the  quantities 
on  the  left  hand  side  of  the  Lagrangian  equation  are  known. 
The  unknown  generalized  force  is  found  by  plugging  in 
suitable  values  into  the  Lagrangian  equation.  For  the 
forward  one.  we  are  given  the  initial  positions,  velocities 
and  generalized  forces.  By  certain  techniques,  which  will 
be  described  later,  we  can  extract  the  unknown  acceleration 
terms  and  solve  for  them.  The  accelerations  obtained  are 
then  integrated  to  find  new  values  of  positions  and 
velocities.       This    step-wise   procedure    is    continued   until 


the    termination    time    is    reached. 

Related  work  in  this  area  has  been  done  in  the 
derivation  of  Lagrange's  equations  of  motion  using  MACSYMA 
[12].  The  objective  of  the  work  was  to  find  the  dynamic 
equations  for  any  general  robot  and  expand  them  into 
symbolic  form.  The  method  used  to  describe  the  system 
kinematics  was  the  Denavi t-Har tenberg  convention  [13]  with 
homogeneous      coordinate      transformations.  The     dynamic 

equations  are  set  up  following  Paul's  derivation  [14] 
which  is  basically  a  series  of  matrix  and  vector 
multiplications.  The  matrix  form  of  the  dynamic  equations 
was  then  expanded  symbolically  using  simple  commands  for 
matrix  multiplication  provided  by  MACSYMA.  To  simplify  the 
expanded  expressions,  it  was  suggested  that  one  use  RATSIMP 
for  algebraic  simplification,  and  TRIGREDUCE  for 
trigonometric  simplification;  one  can  also  define  one's  own 
simplifying  functions  by  using  the  RATSUBST  feature  of 
MACSYMA.  The  benefits  of  this  approach  are  elimination  of 
the  time-consuming  and  error-prone  manual  derivation 
process  and  the  generation  of  equations  which  are  both 
insightful  and  more  computationally  efficient.  The  force 
components  are  derived  individually  as  inertial,  Coriolis, 
and  gravitational  force  terms.  However,  only  the  force 
terms  are  derived  as  explicit  functions  of  other 
generalized     quantities.       [12]      solves     only     the     inverse 


dynamic  problem. 

The  research  in  this  thesis  develops  a  generalized 
technique  for  generating  simulators  for  serial  open  chain 
robots  using  the  symbolic  algebraic  manipulation  language 
REDUCE-3.  Though  the  generality  for  different  robots  is 
obtained  through  symbolic  computing,  it  is  also  known  that 
the  language  is  very  slow  in  arithmetic.  On  the  other 
hand,  FORTRAN  is  relatively  strong  in  numeric  computing  but 
ill-suited  for  symbolic  computing,  the  combined  use  of  both 
languages  is  used  in  the  implementation.  Because  of  the 
generality  and  flexibility,  no  assumptions,  modifications, 
or  simplifications  are  imposed  on  the  original  Lagrangian 
equation,  and  hence  the  solutions  obtained  from  the 
simulator  are  exact.  Symbolic  languages  tend  to  expand  all 
expressions  and  considerable  simplification  is  often 
required  before  these  expressions  can  be  evaluated;  this  is 
one  of  the  difficulties  encountered  in  [12].  In  general, 
we  can  not  guarantee  that  all  the  possible  simplifying 
functions  are  defined  in  the  simulator  generator;  further, 
the  simplification  process  often  takes  too  much  memory 
space  in  computers.  For  certain  very  complicated  problems, 
it  might  just  fail  for  lack  of  memory.  To  avoid 
excessively  lengthy  expressions,  the  individual  terms  in 
the  Lagrangian  equation  are  treated  separately.  By 
expressing  them  as  combinations  of  other  basic  quantities, 
we  can  derive  only  these  basic  expressions  symbolically  and 


set  up  FORTRAN  loops  to  combine  them.  Once  the  lengthy 
expressions  are  avoided,  the  generated  FORTRAN  simulators 
are  ready  for  compilation  and  execution,  and  no  further 
simplification  is  necessary.  The  sensitivity  equations  are 
also  derived  so  that  the  simulators  may  later  be  extended 
to  include  design  optimization  capabilities. 

The  derivation  of  Lagrange's  equations  for  the  forward 
and  inverse  dynamic  problems  are  presented  in  Chapter  2  of 
this  thesis.  The  design  sensitivity  equations,  which  are 
derived  from  the  equations  of  motion,  are  developed  in 
Chapter  3;  the  potential  use  of  these  equations  for 
optimizing  the  design  of  new  robots  and  the  utilization  of 
existing  ones  is  also  discussed.  Chapter  4  describes  the 
implementation  of  the  simulator  generator  with  emphasis  on 
the  interaction  between  the  symbolic  and  numerical  parts  of 
the  code  generation/execution.  In  order  to  test  the 
efficiency  of  the  simulator  generator,  several  examples 
were  run  for  both  the  forward  and  inverse  problems.  some 
of  these  examples,  along  with  the  verification  tests  that 
were  made  to  establish  the  correctness  of  the  results,  are 
presented  in  Chapter  5.  Finally,  the  conclusions  drawn 
from  the  present  work  and  some  recommendations  for  future 
work  in  this  area  are  discussed  in  Chapter  6. 


CHAPTER  2 
SYSTEM  EQUATIONS 

To  analyze  the  response  of  a  dynamic  system,  it  is 
necessary  to  formulate  the  equations  of  motion  for  the 
system  of  interest.  The  equations  of  motion  for  a 
mechanical  system  can  be  obtained  either  from  Newtonian 
mechanics  [15]  or  from  Lagrangian  mechanics  [12,16,17].  The 
Lagrangian  formulation  is  more  popular  because  it  gives  a 
straight  forward  way  of  treating  constrained  dynamic 
systems;  and  most  mechanical  systems  can  be  modelled  in 
terms  of  rigid  bodies  and  constraints. 

In  this  research,  robots  with  a  serial  combination  of 
rigid  links  connected  by  joints  are  of  interest;  such  a 
robot  would  be  completely  described  if  the  links  and  joints 
are  completely  specified.  In  order  to  specify  the  relative 
positions  of  the  links,  each  link  is  associated  with  a 
cartesian  coordinate  that  is  fixed  to  it;  the  homogeneous 
transformation  matrices  between  successive  coordinates  can 
be  found.   These  matrix  transformations  are  of  the  form 

Pi-l=Ai.iPi  (2-1) 

where 

Aj_2  is  a  4x4  transformation  matrix, 

Pj  is  the  homogenious  representation  of  a  vector,  p,  in 
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the    i-^^   coordinate    system. 

The  function  of  a  joint  between  two  links  is  to  provide 
relative  motion  and  is  thus  associated  with  the  degree-of- 
freedom  of  the  system.  Due  to  mechanical  simplicity,  the 
construction  of  most  practical  robots  is  designed  in  such  a 
way  that  there  are  only  single  degree-of -f reedom  joints 
between  the  links.  Two  examples  of  single  degree-of- 
freedom  links  are  the  revolute  joint  and  the  translat ional 
joint.  The  relative  motion  associated  with  revolute  joints 
is  an  angular  rotation  while  that  for  trans lat ional  joints 
is  a  linear  displacement.  By  introducing  suitable 
generalized  coordinates  in  the  Lagrangian  formulation,  we 
can  treat  both  joint  types  in  a  unified  way.  We  can 
construct  a  system  generalized  coordinate  vector,  q, 
consisting  of  the  joint  coordinates;  also,  by 
differentiating  it  with  respect  to  time,  we  can  define  the 
generalized  velocity  vector,  q,  and  the  generalized 
acceleration    vector,     q. 

To  describe  the  links  fully  from  the  point  of  view  of 
dynamics,  it  is  necessary  to  specify  the  position  of  the 
center  of  mass  of  each  link;  and  it  is  convenient  to  give 
this  specification  in  the  corresponding  local  coordinate 
system  since  this  specification  is  invariant  with  respect 
to  motion  of  the  system.  We  also  need  the  mass  of  each  body 
and  the  moment  of  inertia  of  each  body  about  the  local 
coordinate    axis. 


Once      the     above     Information     about      the     system      is 

provided,     it    is    possible    to    express    the   Lagrangian   of    the 

system    in    terms    of     these    parameters.        We    are    presently 

interested      in     both     the     forward     and     inverse     dynamics; 

therefore,    both    sets   of    system   equations   are    required    for 

the    purpose    of    analysis.       In    deriving    the    system    equations. 

we      must     also     ensure      that      they     can     be     conveniently 

implemented    on    the    computer. 

Lagrange's    equation    of    motion    can    be    written   as 

d      dL  dL 

—  {  — ) =   Qj.  j  =  1.2 n  (2-2) 

dt    dqj  dqj 

where 

q -•    is    the    generalized    coordinate, 

q-    is    the    generalized    velocity, 

Q-    is    the    generalized    force, 

n    is    the    degree-of-f reedom    of    the    system, 
and    L    is    the    system   Lagrangian    defined    by 

L=T-V 
where 

T  is  the  system  kinetic  energy, 
and  V  is  the  system  potential  energy. 

To  construct  the  Lagrangian,  the  system  kinetic  and 
potential  energies  should  be  expressed  in  terms  of  the 
generalized  coordinates,  velocities,  and  link  parameters. 
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We  will  first  consider  the  potential  energy.  Since 
spring  forces  are  included  in  the  generalized  force  term  in 
the  equation,  the  only  contribution  to  the  potential  energy 
is  the  gravitational  potential  energy.  Assume  that  the 
base  coordinate  is  established  in  a  way  that  its  y-axis  is 
along  the  vertical;  therefore,  the  gravity  acts  in  the 
negative  y  direction  and  the  potential  energy  of  any  link 
is  given  by 

Vi=migYi  (2-3) 

where 

Yj  is  the  y-coordinate  of  the  center  of  mass  of  link, 

i,   as  measured  from  the  base  coordinate, 

mj  is  the  mass  of  link,  i, 
and  g  is  the  gravitational  acceleration. 

Since  the  coordinates  of  the  center  of  mass  are  given 
with  respect  to  the  local  coordinate,  and  the 
transformation  matrix  transforms  between  two  successive 
frames,  we  can  then  use  a  series  of  transformation  matrices 
to  find  the  coordinate  as  viewed  from  the  base,  given  by 

[Xi,Y..Zi]T  =  AiA2A3...Ai[Xi,yi.Zi]T  (2-4) 

where  the  terms  within  brackets  on  the  left  and  right  hand 
side  are  the  homogeneous  coordinates  of  the  center  of  mass 
in  the  global  and  local  coordinate  systems  respectively. 
The  A*s  are  4x4  joint  transformation  matrices;  their 
products  can  be  written  in  a  simpler  way,  by  defining 
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Tj=AiA2A3...Aj 


(2-5) 


JJi 


Tj   is   then   the   transformation  matrix  from   i-**-"-  local 

coordinate  to  the  base.   The  system  potential  energy  would 

be 

V=  I  Vi=  2  m^gYi  (2-6) 

i      1 

The  motion  of  any  link  can  have  both  trans lat ional  and 

rotational  components;   therefore,   the  kinetic  energy  term 

is  comprised  of  their  respective  contributions.   It  can  be 

writ  ten  as 

T=    —I    [mi(v2^-Hvf    +vfjMIixWfx+IiyWiy+Iiz*iz)] 
2       i  '  ^       ^ 

(2-7) 

where 

v    and    W    are    t rans la t i onal    and    rotational    velocities 

respectively . 

subscripts    x,     y,     and    z    denote    the    component    along    each 

axis  , 

i     is    the    link    number, 
and    the    I's    are    the    centroidal    moments    of    inertia. 

The  translational  velocity  can  be  found  by 
differentiating  Equation  2-4  with  respect  to  time;  noting 
that      Tj      is      an     abbreviation      for      A2A2...Aj,       we      have 

t^ix    Viy    Vj^    l]=[Ti][x^i    y^i    z^i     1]  (2-8) 

and    from    equation    4   we    find 

Tj=AjA2...Aj+AjA2...Aj  +  ...+AjA2...Aj  (2-9) 
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For  the  rotational  velocity  terms,  given  the  moments  of 
inertia  in  the  local  coordinates,  it  is  desirable  to  find 
the  absolute  angular  velocities  for  each  link  expressed  in 
the  local  coordinate  system  of  that  link.  It  is  clear  that 
the  velocity  of  a  link  is  not  kinemat ical ly  dependent  on 
the  links  that  follow  it,  but  only  on  preceding  links  and 
itself;  therefore,  it  involves  only  the  transformation  of 
angular  velocities  from  preceding  links  to  the  local  frame 
of  interest.  We  compute  the  total  angular  velocity  of  link 
,i  .  expressed  in  the  i -^"-  local  coordinate  system 
recursively  as  follows: 

Suppose  the  total  angular  velocity  vector  of  link,  i-1,  is 
known  in  its  local  i-1^-^  coordinate  system,   it  can  be 
converted  into  the  i^-^  coordinate  system  by  premul  t  iply  ing 
this  vector  by  a  transformation  matrix,  given  by 

«i=Rr^*i-i)  (2-10) 

where 

Rj  contains  only  the  rotational  terms  of  Aj. 

w^  is  the  contribution  to  angular  velocity  of  link,  i, 

from  the  preceding  links, 
and  Wj_j  is  the  total  angular  velocity  of  link,  i-1,  given 
by 

Wi_i=Wi_i+Wi_i  (2-11) 

where 

Wj_j  is  the  angular  velocity  of  link  .  i-1,  itself. 


13 


—  1  T 

Since  Rj  Is  an  orthogonal  matrix.  Rj   is  just  Rj  which 

is  easy  to  construct  symbolically.   Thus  for  any  link,  i. 

we  can  find  the  total  angular  velocity  in  the  following 

way.   Starting  from  the  first  link,  we  find  the  total 

angular  velocity  of  link,  1,  expressed  in  local  coordinate 

of  link,  2.  This  is  then  added  to  the  angular  velocity  due 

to  relative  motion  between  link,  1  and  2,  to  obtain  the 

total  angular  velocity  of  link,  2.   We  proceed  in  the  same 

way  until  the  total  angular  velocity  of  link,  i,  is  found. 

Note  that  for  a  t rans la t i onal  joint,  the  relative 
angular  velocity  between  the  connected  bodies  is  zero. 
Once  the  trans lat ional  and  rotational  velocities  are  known, 
it  is  easy  to  explicitly  derive  the  system  Lagrangian  in 
symbolic  form  by  using  simple  loops  in  a  REDUCE  program. 
However,  this  is  not  advisable  because  the  velocity  terms 
may  be  quite  complicated  and  squaring  them  could  lead  to 
expressions  that  are  too  large  to  handle.  What  we  need 
presently  are  not  the  expressions  for  kinetic  energy  or  the 
Lagrangian;  rather,  we  are  interested  in  the  derivatives  of 
the  Lagrangian  with  respect  to  generalized  coordinates  and 
velocities.  We  must  therefore  develop  formulas  for  these 
quantities  which  do  not  require  the  evaluation  of 
unreasonably  large  expressions. 

The  Lagrangian  equation  can  be  expanded  term  by  term 
in  the  following  way: 
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dL  d{T-V) 

dqj  aqj 

dv.  aw.  dYj 

=  Z   miVi    — -    +  I    liWj    — -    -  Z    m^g    —  (2-12) 

i  dqj  i  &qj  i  3qj 


and 


d      dL  d      a{T-V)         d      hT 

__(___)=    ____    =    __(„__)  (2-13) 

dt    9q^         dt       3qj  dt    3qj 

The    potential     energy     term     is    dropped     from     the    above 

equation    because    it    is    independent    of    velocities. 

d      dT                 1      d      d  Z    (m.vf+I.W?) 
__(___)    = |;__i ] 

dt    dq-  2       dt  flq. 

dv.    dv.  d      dv.  ^  dW.    dW.  d      SW. 

=  I   „,.,;__i    ..i    +    vj    -(--)]    +  I    liC "    +    Wi    -(---)] 

1  dt       dq^  dt    dq  •  i  dt       9q •  dt    dq^ 

(2-14) 

Now  we  can  find  the  Lagrangian  equation  by  subtracting 
Equation  2-12  from  Equation  2-14  and  setting  it  equal  to 
the  generalized  force.  The  Lagrangian  equation  is  then  in 
a    form     that    can    be    handled    by    REDUCE. 

For  the  inverse  dynamic  problem,  we  want  to  find  the 
generalized  forces  required  to  satisfy  prescribed 
requirements  on  the  displacements,  velocities,  and 
accelerations.  We  see  that  the  unknown  generalized  force 
appears  on  the  right  hand  side  of  Equation  2-2,  and  that 
all  the  terms  on  the  left  are  known  quantities;  thus  the 
evaluation    of     the    forces    is    straight    forward. 
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For  the  forward  dynamic  problem,  we  wish  to  calculate 
the  displacements,  velocities,  and  accelerations  given  the 
generalized  forces  and  the  initial  displacements  and 
velocities.  The     unknown     accelerations     are      solved 

algebraically  from  the  equations  of  motion,  and  then 
integrated  to  obtain  velocities  and  positions.  However, 
the    accelerations    for    which    we    must    solve   are    mixed    with 

dVi 

other  known   quantities   on   the  left  side  in  the and 

dt 


dWi 


dt 


terms.    Again,  there  is  a   perfect   similarity  between 


trans lat ional  and  rotational  accelerations,  because  of  the 
use  of  the  generalized  coordinates.  Since  the  accelerations 
are  obtained  from  multiplication  of  three  matrices  and 
there  is  only  a  first  power  of  q.  it  follows  that  the 
multiples  are  linear  combination  of  accelerations.  All  the 
terms  containing  accelerations  are  then  extracted,  while 
the  other  known  quantities  are  moved  to  the  right  hand 
side.  Once  this  is  done,  the  resulting  set  of  linear 
equations  can  be  solved  for  the  accelerations.  Finally, 
the  solved  accelerations  are  integrated  in  a  given  time 
interval  to  find  the  velocity  and  the  displacement  for  the 
next  step  by  using  the  Runge-Kutta  integration  scheme. 
Provided  the  generalized  forces  for  the  next  step,  the 
process  is  repeated  until   the  end  of  the  simulation. 
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CHAPTER  3 
DESIGN  SENSITIVITY  ANALYSIS 

The  utility  of  a  simulator  in  the  design  process  is 
greatly  enhanced  if  the  simulator  has  the  ability  to 
automatically  optimize  the  system  design.  It  is  envisioned 
that  the  methods  developed  in  the  preceding  chapter  for 
analysis  will  be  extended  to  include  optimization 
capabilities.  This  will  enable  the  use  of  the  simulator 
for  optimizing  the  design  of  new  robots  as  well  optimizing 
task-planning  for  existing  ones. 

Symbolic  computing  is  particularly  useful  in  an 
optimization  environment  because  optimization  problems  are 
generally  more  difficult  to  standardize;  there  can  be  a 
wide  variety  of  choices  of  design  parameters  and  optimality 
criteria.  Thus,  it  is  desirable  to  maintain  the  symbolic 
computing  approach  while  extending  our  analysis  method  into 
the  realm  of  optimization. 

A  method  of  performing  this  extension  in  a  manner  that 
will  permit  convenient  implementation  in  a  symbolic 
manipulation  language  will  be  developed  in  this  chapter. 

For  the  purpose  of  optimizing  a  design,  the  design 
problem  should  be  expressed  in  a  standard  format.  The 
standard  problem  is  stated  in  the  following  manner. 
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Suppose  the  design  of  the  system  is  specified  by  a 
vector,  b,  of  r  design  variables,  i.e. 

b=[bi.b2 bj.]"^  (3-1) 

Find  thedesign,b,  which  minimizes  a  specified  cost 
function,     Fq,     subject    to    constraints 

Fj=0,  1  =  1.2 m 


(3-2) 


Fj^O.  i  =  m+l.m  +  2 n 


The  optimal  design  process  solves  this  standard  problem 
by  iterative  improvement  of  the  design  vector,  b.  Each 
iteration  of  this  solution  process  consists  of  three 
di  s  t  inc  t  s  teps : 

1.  System  analysis:  The  behavior  of  the  system  is  analyzed 
at  the  current  design.  The  cost  and  constraint 
function  values  are  also  computed. 

2.  Design  sensitivity  analysis:  The  derivatives  of  the 
cost  and  constraint  functions  with  respect  to  the 
design  variables  are  calculated  in  this  step. 

3.  Optimization/design  update:  The  cost  and  constraint 
function  values  from  the  system  analysis  as  well  as  the 
derivatives  obtained  from  design  sensitivity  are 
supplied  to  an  optimizing  algorithm.  The  algorithm 
computes  the  required  change  in  design  and  proceeds  to 
the  next  iteration  unless  a  predefined  convergence 
criterion  is  satisfied. 

The  system  analysis  is  done  by  solving  the  system 
equations  developed  in  Chapter  2.   The  design  sensitivity 
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analysis      is     discussed     in     this     chapter     for     the     future 
purpose    of    optimizing    and    updating    the    design. 

The  system  equations  of  motion  derived  in  the  preceding 
chapter  are  dependent  on  the  design  vector,  b.  of  the 
system.  The  optimizing  process  is  to  choose  suitable 
values  for  the  components  of  the  design  vector  in  order  to 
minimize  a  cost  function,  subject  to  the  given  constraints. 
We  assume  that  the  cost  and  constraint  functions.  F^, 
depend  on  time  (t),  position  (q),  velocity  (q), 
acceleration  (q),  joint  generalized  force(Q),  and 
design(b).       That    is,     they    are    of    the    form 

Fi=Fi(t.q,4.q.Q.b),  i=0,l n  (3-3) 

Most  optimization  algorithms  require  the  derivatives  of 
cost  and  constraint  functions  with  respect  to  design.  The 
aim  of  first  order  design  sensitivity  analysis  is  to 
evaluate  these  derivatives  to  first  order.  These 
derivatives  relate  variations  in  cost/constraints  function 
values  to  those  in  the  design  parameters  through  the 
equation : 

Sf  =  I'^Sb  (3-4) 

where   1   is   the  matrix  (nxr)   of  design  sensitivity 
coef  f  icients . 

In  order  to  find  the  design  sensitivity  coefficients, 
we  have  to  take  the  variations  of  the  F^'s.  which  are  of 
the  form 
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dF      dF      dF      dF      dF 
St    =  --Sb  +  --6q  +  --&q  +  --S*q  +  --8Q  (3-5) 

db      dq      dq      dq      dQ 

Since  state,   velocity,  acceleration  and  generalized  force 

are  all  implicit  functions  of  the  design,  it  follows  that 

dq  dq  dq  dQ 

dq  =  --Sb.   dq  =  --db.   Sq  =  --db.   dQ  =  --8b      (3-6) 
db  db  db  db 

Equation  3-6  can  be  substituted  into  Equation  3-5  to  obtain 

dF      dF  dq      dF  dq 
dF=  — db+  --  --db  +  --  --db 
db      dq  db      dq  db 

dF  dq      dF  dQ 
+  —  —6b  +  --  --6b  (3-7) 

dq  db      dQ  db 

comparing  Equations  3-7  and  3-4,  we  conclude  that 

dF    dF  dq     dF  dq     dF  dq    dF  dQ 
iT  =  -_  +  __  __  +____+  __  __  +  __  __  (3-8) 

db    dq  db     dq  db     dq  db    dQ  db 

As  with  the  dynamics  analysis,  the  sensitivity  analysis 
also  has  both  forward  and  inverse  cases.  For  the  inverse 
sensitivity  analysis,  the  only  unknown  on  the  right  hand 

aq 

side  of  Equation  3-8  is  which  can  be  found  by 

db 

differentiating  the  Lagrange  equation  with  respect  to  the 

design,  b.  For  the  forward  dynamics  problem,  we  must  have 

dq     dq   dQ         dQ 

ini  t  ia  1  va  lues  of , . .  and at  each  step. 

db     db   db         db 

Then.   by  taking  the  first  variation  of  the  Lagrangian 

equation,   we   can   obtain   a   second   order   ordinary 
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differential     equations     in    .      which    can    be    numerically 

db 

integrated  along  with  the  equations  of  motion.   At  any 

dq    dq        dq 

instant,  the  values  of . ,  and thus  obtained  can 

db    db       dh 

be  substituted  into  Equation  3-8  to  obtain  the  design 

sensi  t ivi  ty . 

dQ 

We  thus  see  that  the   unknown   terms   are, ,  for  the 

db 
dq     dq         dq 

inverseanaly s i s ,  and  ,  ,  and  for  the  forward 

db     db        db 

analysis.   It  is  desirable  to  derive  these  terms  so  that 

they  are  ready  for  implementation  in  a  computer  code.   When 

differentiating  a  series  of  matrices,  followed  by  a  vector, 

with  respect  to  another  vector,  it  is  possible  to  avoid 

three  dimensional  tensors  within  the  process  by  using  the 

mathematical  technique  described  below. 

Suppose  we  want  to  find 

i—  (3-9) 

db 

where 

Aj  is  a  two  dimensional  matrix, 

and  b  and  x  are  vectors. 

The  technique  is  given  by 

d{Aix)      d(A,x)  dx 

i...    = i___  +   Aj (3-10) 

db  db  db 
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The  tilde  (~)  in  Equation  3-10  means  that  the  term 
beneath  it  is  treated  as  constant  when  doing 
differentiation,  but  after  the  process,  it  is  again 
recognized    as    its    original    nature. 

For  a  more  complicated  case,  we  can  generalize  the 
procedure    as    follows. 


--—  (A1A2A3.  .  .A.x) 


(3-11) 


(Nest)i^  =  Dj^  {Nest)i^+j  + 


can  be  expressed  into  a  form  of  a  nest  within  another  by 

^{Dk_Bk)_ 
db 

k=j.j-l 1  (3-12) 

where 

Dj^  is  the  dummy  that  represents  transformation  matrices 

A's. 
and  Bj^  is  the  multiple  from  the  k+1^^  matrix  to  the  end  in 
(A1A2A3. . .AjX). 

The  iteration  starts  with  k=j  and  ends  with  k=l  and 

(Nest)^  is  the  result  of  the  derivation.   To  start  the 

iteration.  (Nest)^^j^  must  be  given.  Therefore,  we  define 

dx 

(3-13) 


(Nest)j^.i  = 


db 


Equations  3-12  and  3-13  altogether  is  the  solution  of 


Equation  3-11.   Base  on  this,  we  can  find  the  terms 


db 
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and  that  are  necessary  for  the  forward  and  inverse 

db 

sensitivity  analyses.   The  key  to  both  terms,  as  stated 

previously,  is  to  start  with  taking  the  derivative  of  the 

Lagrange  equation  with  respect  to  the  design  variable. 

d   d   dL       aL 

-[--{-:-)  -  (  — )] 

3b  dt  dqj      ftqj 

3m,.  dv^  dv,-    _     d  dv  ^    dv ..         dv .-  d  dv  ^ 

=  2  -'-    -'-    -'-    +  Z  m,  -  {  -^  )  -^  +  £  mj  -^  -  (  --'-    ) 
j  db   dt   ^qj    j     db    dt     dq^  j     dt   db    dq^ 

^  dm..    dY-.  d     dY. 

+  Z  — ^  g  — ^  +  Z  THj  g  —  (  — ^  )  (3-14) 

j  db     dqj    j       db    dq^ 

For  the  inverse  sensitivity  problem.  Equation  3-14  is 

dQ 

just  ,  the  unknown.   To  evaluate  3-14,  some  of  the  terms 

db 

d   dvj. 
are   seen   to   be   very   complicated,   such  as   — ( ) 

db  dt 

d   dVj. 

and  --( ),   from  the  experience  in  dealing  with  the 

db  d^i 

dynamics  analysis.   The  other  terms  are  either  quite  simple 

for  REDUCE  to  handle  or  has  been  done  in  the  dynamic 

problem . 

Following   the  convention  described   in   the  second 

dv  ^ 

chapter,  we  can  expand  in  the  following  manner. 

dt 
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3   dv.     d    ..  .  . 

__(__£)  =  --[(AjAg. . .Ajr+AjAg. . .Ajr+. . .+A1A2- •  AjT) 

9b  dt      dh 

•  •  .  ■  ■•     » 

+(A2A2- • .A.r+AjA2. . •Ajr+. . .+AjA2. •  Ajf) 


+  (A2A2-  •  .A.r+AjA2-  •  .A  ■r+.  .  .+A2A2.  •  -Ar)] 

(3-15) 

Let's  take  every  AjA2---A^r  as  a  group.   The  previous 

described  technique  can  be  applied  to  obtain  the  derivative 

for  each  group.   Then,  by  summing  the  derivatives  for  all 

the  n  X  n  groups,  we  can  find  the  solution. 

dv 


We  can  also  expand  similarly 


d   dv.     d         dk. 

__(__J)  =  __{  ___  A2...Ajr 


db    bq^  dh      dq^ 


dA2 

+Ai  ---  ...Ajr 


dqi 


aA. 

+A,A9. . .  — -  r)      (3-16) 
dqi 

Since  q^  will  exist  in  A-  only  when  i=j,  some  of  the 

derivatives  in  Equation  3-16  are  always  zero  and  hence  not 

shown.   Actually,  the  only  non  zero  term  in  Equation  3-16 

is  the  i^^  group.   Again,  the  derivative  can  be  evaluated 

in  the  previous  way. 


24 


Once  these  complicated  terms  are  obtained,  along  with 

simpler  terms  derived  by  REDUCE,  we  are  ready  to  calculate 

5Q 

using  Equation  3-14.   It  is  then  substituted  back  into 

db 

8F 

Equation  3-8  to  obtain  which  is  the  objective  of  the 

5b 

inverse  sensitivity  analysis. 

For  the  forward  sensitivity  analysis,  we  also  have  to 

follow  Equation  3-14.   This  time,   we  are  interested  in 

dq 

isolating  the  unknown .   From  Equation  3-14.   it  is 

dh 

dq                           9   dv  • 
apparent  that is  mixed  with  others  in  — ( )  only,  and 

db  db  dt 

hence  we  are  to  concentrate  on  this  term. 

The  term,  ,  can  be  expanded  as  in  Equation  3-15.  The 

dt 

q's     exist     only     in     the    diagonal     groups.        Therefore,      the 

diagonal    groups    shoud    be    separated    into    two    parts.       One 

with    q    and    the    other    without,    as    follows. 

[Ai]=[AA]q.+[AB]  (3-17) 

Accordingly,      Equation    3-15    becomes     the     combination    of 

two ,  i.e. 

d     ..g  .  .  .       . 

—  [(AjA2. . .A.r+AjA2. • •Ajr+. . .+A2A2. • -A-r) 
^b    .  .  .  .     . 

+(AjA2. . .AjT+AjAg. . •Ajr+. . .+AjA2. •  Ar) 


+(AiA2. . .Ajr+AiA2. • •Ajr+. . .+A1A2. .  A^r)]       (3-18) 
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and 


ob 


.Ajr 


+AiAA...Ajr 


+AiA2...A5r]  (3-19) 

Equation  3-18  is  composed  of  constants  and  is  to  be 

moved  to  the  right  hand  side  of  the  Lagrange  equation.   The 

i^^  group  in  Equation  3-19  corresponds  to  the  contribution 


d      dv 


J, 


of   —  { )   to   the   coefficient   of   q^.    i  =  l,2 j. 

db  dt 

Equations  3-18  and  3-19  are  evaluated  by  the  previous 

method  again. 

d      dvj 

After  we  go  through  every  --( ),  j  =  l,2 n,  we  can 

db  dt 

dqj 

set  up  a  coefficient  matrix  for ,  while  moving  and 

db 

combining  all  the  known  quantities  in  Equation  3-14  with 

the  generalized  force  on  the  right  side  of  the  Lagrange 

equation.    We  can   solve  the  n  simultaneous   equations  for 

,  1  =  1,2 n.   They  are  then  substituted  into  Equation 

db 

5f               dq 
3-8  to  obtain  the  sensitivity, .  Note  that are  the 

5b  db 

only  unknowns  in  that  equation. 
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dq       dq 

Finally.  and  would  be  predicted  by  integrating 

db      db 

dq 
at  this  time,  with  updated  values  of  generalized  forces 

db 

and    the    state,    velocity,    and    acceleration    obtained    from 
dynamic    analysis.       We    are    able    to    resume    the    iteration    for 

dq 

another  until  the  end  of  the  simulation. 

db 
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CHAPTER    4 
IMPLEMENTATION 

In    the    implementation    of    system    equations,     the    symbolic 
manipulation    language    was    chosen    for    following    benefits. 

1.  All  the  derivatives  involved  in  the  development  of  the 
system  equation  are  found  automatically  and  then 
calculated  by  the  preprocessor  so  that  no  further  user 
effort  is  necessary.  On  the  other  hand,  if  an  all- 
FORTRAN  program  is  used,  the  user  must  supply 
subroutines    for    evaluating    the    required    values. 

2.  Most  of  the  errors  that  occur  in  the  development  and 
use  of  software  arise  from  mistakes  made  by  either  the 
programmer  or  the  user.  By  using  computers  to 
symbolically  derive  the  equations,  the  chance  of 
programmer  error  is  greatly  reduced.  Furthermore,  by 
minimizing  the  input  required  from  the  user,  the 
probability  of  user  error  is  also  lowered.  These  two 
factors    make    REDUCE   based    programs    more    reliable. 

3.  The  FORTRAN  subroutines  generated  by  REDUCE  are 
p r ob 1 e m - d ep end en t  because  they  are  constructed 
specifically  for  the  problem  that  is  being  solved. 
They  are  much  more  efficient  than  the  general  purpose 
subroutines.       Therefore,     the    regular     trade-off    between 
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the  generality  and  its  efficiency  does  not  apply  to 

REDUCE  based  software. 
4.   Explicit  formulation  of  expressions  can  be  obtained  in 

FORTRAN  code.   These  could  be  of  interest  to  the 

analysts  or  designers  and  of  help  in  understanding  the 

behavior  of  the  system. 

The  system  equations  of  motion  have  been  derived  for 
both  the  forward  and  the  inverse  dynamic  problem  in  Chapter 
2.  All  the  terms  in  these  equations  can  be  found  by  using 
simple  differentiation,  and  the  four  basic  arithmetic 
operations  in  REDUCE.  The  basic  solution  procedure  to  be 
followed  is  to  feed  a  system  description  to  the  REDUCE 
program  and  have  the  program  generate  a  FORTRAN  simulator 
tailored  to  the  specific  robot  of  interest.  We  can  then 
feed  this  FORTRAN  program  with  a  set  of  numeric  data  and 
solve  the  equations  of  motion  numerically.  The  two  step 
approach  was  used  successfully  in  the  automatic  generation 
of  optimizing  simulators  for  general  constrained  planar 
mechanical  systems  [18].  The  rationale  for  this  approach 
is  that  a  symbolic  manipulation  language  has  a  much  higher 
degree  of  generality  and  flexibility,  but  at  the  same  time, 
is  very  slow  in  arithmetic.  To  remedy  this  handicap,  many 
symbolic  languages,  including  REDUCE,  provide  a  facility 
for  writing  out  the  expressions  they  generate  in  the  form 
of  FORTRAN  statements.  This  approach  provides  broad 
generality  and  flexibility  in  the  symbolic  computing  step 
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while    the    resulting    FORTRAN   program    reflects    the    efficiency 

and    problem    specificity    of    a    special    purpose    simulator. 

The    number    of    additions    and    multiplications    associated 

with     the    Lagrange    dynamic    formulation    exhibits    quartic 

variation    with    respect    to    the    degrees    of    freedom    of    the 

system.        If    we    write    the    system    Equations    explicitly    as 

symbolic     languages     usually     do.      one     can     imagine     how 

cumbersome     the     expressions     would    be     for     a     system     with 

several     degrees     of     freedom.        Further,     by    expanding    all 

terms,    we   are   placing    too    many    computational    demands    on    the 

resulting    FORTRAN    program     since     the    expansions    are     rarely 

computationally     efficient.         To     avoid     both     the     huge 

expressions    and    the    redundant    evaluations,     we    can    write 

out    only    the    frequently  used   basic    terms.       Evaluate    them 

once  and  keep  them   in  memory,    then  we  can  use  FORTRAN  loops 

to    do     the    multiplications    and    additions    on    these    known 

quantities.       In    the    progress    of    this    research,     it    was    found 

that    most    of    these    terms    are    quite    short   and   can    reasonably 

be    handled    by    the    FORTRAN    program    with    every    term    expanded. 

dv    3v         d   dv 

Only  certain  terms  such  as ,   ,   and  --( )  ,need 

dt    dq         dt  dq 

special  treatment. 

It  can  be  seen  that  the  velocity  plays  an  important 

role  in  complicating  the  problem;   therefore,   it  is 

desirable  to  express  the  velocity  term  in  a  way  that  can  be 
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handled  by  FORTRAN  loops.  The  velocity  of  the  center  of 
mass  of  link  i  is  obtained  by  differentiation  of  the 
position  of  the  center  of  mass  with  respect  to  time,  i.e. 

Vi=Ri     .   i  =  1.2 n  (4-1) 

Note  that  Vj    and  Rj  are  measured  in  the  base  frame. 

The  position  of  the  center  of  mass  is  given  in  Equation  2-4 

by 

Ri=AiA2...Ai[ri]'r  (4-2) 

The  velocity  of  the  center  of  mass  is  then 

d 
Vi=— {AiA2...Ai[ri]T^}  (4-3) 

dt 

Note  that  [r^]  is  constant  in  the  local  frame 

Vi=[AiA2. . .Ai+AjAg. . .Ai+. ..+A1A2. . .AjJCrj] 


(4-4) 


now  we  want  to  find  the  time  derivative  of  Vj. 

dvj      ..  .  .        .       . 

=  [(A2A2. . . AJ+A2A2-  • AJ+A2A2. • • Aj) 

dt 

+(AiA2. . .Ai+AiA2. - .A^+A^k2. ■ -Aj) 


+(AiA2. . .A^+AjAg. . .AJ+A1A2. - •Ai)][ri] 

.1  =  1.2 n  (4-5) 

The  Equation  4-5  is  now  in  a  form  that  is  suitable  for 
programming  in  FORTRAN  loops.  The  key  is  to  identify  which 
matrices  should  be  used  in  the  different  matrix  products. 
If  we  view  the  terms  in  Equation  4-5  as  consisting  of 
distinct  groups,   it  is  clear  that  it  is  made  up  of  ixi 
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groups,  while  each  group  consists  of  1  transformation 
matrices.  We  can  Index  the  group  by  means  of  two  Indices, 
namely  1  and  m.  and  the  term  number  within  each  group  by  k. 
Once  each  term  has  been  Indexed  In  this  manner,  the 
following  logic  holds: 
when   l=m 

if    k=l     then    we    take    Aj^, 

for  other  k,  we  take  Aj^. 
when  l^m 

if  k=l  or  m  then  we  take  A-^. 

for  other  k,  we  take  k^. 

Similarly  for  ,  we  can  write 

dv^       dAj  dA2 

---  =   (  ---  Aa-.-A^r  +  Ai  -- 


dq 


dq 


k^T+.. .+A1A2. 


uAi   _  "^2 

+(  Ag.-A^r+Aj  ...Ajr+ 

dqj  dqj 


+  A1A2 


SA, 
r) 

r) 

aqj. 


dAi        .        ^^2      •  ^^i 

+(  A2- •  A^r+Aj  . . . A^r+. . . +AjA2- • •  r) 

aqj  dqj.  9qj. 

(4-6) 

Since  Aj  contains  only  the  local  generalized  coordinate 

q^,  we  can  further  simplify  the  Equation  by  setting  the 

3Aj^     dk^ 

or terms  equal  to  zero  for  kiij.   It  is  easily  seen 

dqj     aqj. 
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that  only  the  j-"^  column  of  groups  are  non-zero  for 

Using  the  same  indexing  scheme  for  1,  m,  and  k  as  before, 
the  logic  for  evaluating  this  term  is  as  follows: 
For  any  fixed  j,  only  m=j  is  considered 
when  l=m 

if  k=l  then  we  take 

for   other   k,    we    take   Ai^. 
when  l^m 

if  k  =  m   then  we   take 


dqy 


if    k=l     then    we    take    Ai 


for    other   k,    we    take   Ai 


d   9vj 

--{ )  can  be  expanded  as 

d  t  dq  . 


d   dvj     d    dAj 


__(___)  =  __(  ___  Ag.-.Ajr 

dAo 


dt  dq.     dt   dq . 

J  J 


+Ai  ---  ...Ajr 
dqj. 


6A, 


+A1A2...  — -  r)      (4-7) 

aq, 


Note  that  the  terms  are  dropped  because  they  are 

dqj 
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always  zero.   The  remaining  terms  are  not  zero  only 

^^i  th 

when  k  =  j:  therefore, corresponds  to  the  j-*"*^  group  in 


Equation  4-7. 

A  further  differentiation  with  respect  to  time  yields 

d   dAj             dAj  .  dAj 

[  —  ( )  Ag-.-A^rH- Ag.  .  .  Ajr+.  .  .+  ---  A2..-Ajr] 


dt  dq 
dA. 


8qj 

d   dAg 


OCj 


bAc 


+[Aj  ---  ...A^r+Aj  —(-7-)  . . .Ajr+. . .+Ai  ---  ...A^r] 


dq 


dt  dq  . 


2>q 


dAj      .      dAj  d   ftAj 


+[AjA2...  r+AjAg ...  r+ . . . +A jA2 • . .  — (-7-)r] 


dq 


d   dv. 


dq 


dt  hq 


(4-8) 


Now,  for  --( ),  the  only  non  zero  terms  correspond  to 

dt  8q  . 

the  j-"^  row  of  groups.   With  1,  m,  and  k  defined  as  before, 

we  have  the  following  logic. 

For  any  fixed  j,  j  =  l,2 n,  only  l  =  j  need  be  considered 

when  m  =  l 

d  dk^ 

if  k=m  then  we  take  --( ) 

dt  dctj  . 

for  other  k,  we  take  Aj^. 

when  m^l 


dAj^ 

if  k=l  then  we  take 

dqj. 
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if  k  =  m  then  we  take  Aj^, 

for  other  k,  we  take  Aj^. 

Studying  these  above  expressions,   we  see  that  the 

repeatedly  appearing  basic  terms  are  the  A's,  A's,  A's, 

dAj^        dAj^  d      bk^  bky. 

— ( ),      and        .     These     expressions     are 

dqj  3qj  dt    aqj  dqj 

reasonably    short    and   can  be   explicitly    written   out    in    the 

first      part      of      the     FORTRAN     program.         The      matrix 

multiplication   and    the    logic    for    evaluating    the   appropriate 

matrix    products    are    then    programmed    in    FORTRAN    do    loops 

that     are    also    generated    by     the     symbolic     language.        The 

other    terms    in    the    system    equations    can   also   be    expanded    in 

the    same    way;     however,     it    seems    to    be    unnecessary    since 

many    elements    in    the    matrices    are    zero.        It    would   be    easier 

to    write    them    explicitly    then    to   do    multiplications    on    lots 

of    zeros. 

As  stated  in  previous  chapters,  the  simulator  developed 
can  simulate  both  forward  and  inverse  dynamics  of  robot 
manipulators;  therefore , there  are  two  separate  programs  for 
the  two  different  purposes. 

For  the  inverse  dynamics  problem,  the  FORTRAN  program 
starts  by  reading  in  the  dimensions  and  specifications  of 
the  robot  by  calling  an  input  subroutine.  Next,  we  call  a 
subroutine  to  initialize  the  joint  coordinates,  velocities, 
and  accelerations.   After  this,  suitable  loops  are  set  up 
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to    obtain    the    basic    terms    described    earlier    and    combine 

them    to    obtain    the    required    generalized    forces. 

The    basic     layout    of     the    forward    dynamic    simulation    is 

pretty    close    to    the    inverse   one,    except    that    the   big    loop 

now    goes    into    a    subroutine    and    becomes    a    function.       The 

function    in    turn    is    required    by    DVERK.     a    subroutine    in    the 

IMSL     library,      that     does     integration    using     RUNGE-KUTTA 

method.         Also,       what      we     are     interested     now      is      the 

dv 

accelerations;   therefore,   the  term  developed  for  the 

dt 

inverse  problem  is  no  longer  necessary.    It  is  almost 

impossible  to  extract  the  coefficients  of  the  accelerations 

numerically,  fortunately,  the  REDUCE  provides  powerful 

commands  that  enables  us  to  seperate  the  terms  with  and 

without  accelerations.    The  coefficients  are  further 

extracted  from  the  terms  with  accelerations  to  form  an  nxn 

matrix  with  known  quantities,  while  all  the  terms  without 

accelerations  are  moved  to  the  right  hand  side  to  be 

substracted  from  the  generalized  forces.   These  terms  then 

form  a  vector  with  known  quantities.    The  unknown 

acceleration  vector  is  also  established,   then  the  problem 

becomes  solving  an  n-variable  simultaneous  equations.   This 

is  easily  done  by  calling  LEQT2F,  another  IMSL  subroutine. 

Once   the   accelerations   of   links   subjected   to   the 

corresponding  joint  force  are  solved,   they  are  then 

integrated  within  specific  time  step  to  find  new  velocities 
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and    positions    at    next    step.       The    generalized    forces    are 
updated    before    loop    for    the    next    step. 
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CHAPTER  5 
NUMERICAL  EXAMPLES 

The  techniques  for  forward  and  inverse  dynamics  that 
were  described  in  Chapter  2  were  implemented  in  a  REDUCE 
program  and  this  code  was  used  to  generate  simulators  for 
several  example  problems.  These  simulators  were  then 
compiled  and  run  to  obtain  the  system  response. 

The  results  obtained  for  some  of  these  examples  are 
presented  in  this  chapter.  All  the  examples  were  verified 
by  two  independent  methods: 

1)  Work-energy  balance:  A  running  total  of  the  cumulative 
work  done  on  the  system  and  the  total  energy  of  the 
system  is  computed.  This  information  is  used  to  check 
whether  the  results  obtained  are  consistent  with  the 
work-energy  balance  equation  that  must  hold  for  the 
system.  This  check  is  performed  independently  for  the 
forward  and  inverse  dynamics  problems.  The  trapezoid 
approach  is  used  to  find  the  work  done  on  each  joint. 
Though  it  is  not  the  precise  solution,  just  enough  for 
the  purpose  of  checking  the  results.  In  all  the 
examples  that  were  run,  the  check  was  satisfied  within 
the  limits  of  acceptable  numerical  error. 


38 


2)  Cross-check  between  forward  and  inverse  dynamics:  For 
every  system  that  was  considered,  a  known  force  input 
is  supplied  to  the  forward  dynamics  problem  and  the 
resulting  accelerations,  velocities,  and  positions  are 
obtained.  These  results  are  then  used  as  inputs  to  the 
inverse  dynamics  problem  for  the  same  system  and  the 
forces  required  to  produce  this  motion  are  computed. 
These  can  then  be  checked  against  the  original  forces 
that  were  input  to  the  forward  dynamics  problem. 
Ideally,  we  should  obtain  exactly  the  same  forces  in 
the  two  cases;  in  practice,  it  was  found  that  they 
agreed  quite  accurately,  i.e.  within  the  limits  of 
numerical  error. 
In  all  the  examples  presented  on  the  following  pages, 

the  results  of  the  work-energy  balance  and  the  cross  check 

are  presented  along  with  the  system  response  that  was 

obtained  by  running  the  simulator. 


39 


FXAMPT  F  1  : 3-R  RQRQT 

The  initial  position  or  home  position  and  coordinate 
systems  are  shown  in  figrue  1.  The  only  load  is  the  self- 
weight  of  the  members;  the  simulation  time  is  from  0.00  to 
0.18  second . 

The  input  data  is  given  below,  and  all  units  are  in 
MKS. 

Lengths : 

Li=1.5:   L2=1.0;   L3=0.75: 

Masses : 

mi=45;   m2=30;   m2=25; 

Moments  of Inertia^ 


ly 


l2x=0:   l2y=0 
l3x=0:   l3y=0 


Il^=8.4375: 

Il^=2.5: 

l3^=1.17: 


Centers    nf    mass: 

ri^=0.75:       riy=0:        ri^=0: 


r2x=0-50;        r9„=0: 
r3^=0.375: 


2y 

■3y 


2z 


=0; 


=0; 


3z 


=0: 
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(1-1)  Forward  dynamir  prnblem: 
Initial  rnndit.inns: 

qj=0;   q2=0:   ^3=0: 

4^=0.2094;   q2="013^^:   q3=-0.0698; 
GpnaraliTied  fnrcf  fimrtinn  at  each  instant: 

Qj=  lO.Oxt^^^+e.O; 

Q2=  8.0Ht^/2+5.0; 

Q3=    6.0«t^'^2+4.0: 
where    t    is    the    simulation    time. 


RfLSH 

Irs: 
Fnr  link  1 : 

t 

^1 

^1 

0.00 

O.OOOOe+00 

0.2094 

0.02 

0.4200e-02 

0.2106 

0.04 

0.8424e-02 

0.2118 

0.06 

0.1267e-01 

0.2131 

0.08 

0. 1695e-01 

0.2145 

0.  10 

0.2125e-01 

0.2159 

0.  12 

0.2559e-01 

0.2174 

0.  14 

0.2995e-01 

0.2188 

0.  16 

0.3434e-01 

0.2203 

0.  18 

0.3876e-01 

0.2217 

^1 

0.5179e-01 
0.6123e-01 
0.6494e-01 
0.6764e-01 
0.6974e-01 
0.7135e-01 
0.7249e-01 
0.7313e-01 
0.7318e-01 
0.7254e-01 
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For  link  ?.: 

t 

<12 

^2 

0.00 

O.OOOOe+00 

-0. 

1396 

0.02 

-0.2919e-02 

-0. 

1526 

0.04 

-0.6114e-02 

-0. 

1670 

0.06 

-0.9604e-02 

-0, 

1822 

0.08 

-0. 1341e-01 

-0. 

,  1980 

0.  10 

-0. 1753e-01 

-0. 

,2144 

0.  12 

-0.2198e-01 

-0 

.2313 

0.  14 

-0.2678e-01 

-0 

.2485 

0.  16 

-0.3193e-01 

-0 

.2660 

0.  18 

-0.3742e-01 

-0 

.2838 

^12 
-0.5845 

-0.6994 

-0.7462 

-0.7815 

-0.8104 

-0.8348 

-0.8554 

-0.8724 

-0.8855 

-0.8945 


For  link  3: 

t 

^3 

^3 

'^3 

0.00 

O.OOOOe+00 

-0.6981e- 

-01 

2.2970 

0.02 

-0.8966e-03 

-0.1862e- 

-01 

2.7667 

0.04 

-0.7056e-03 

0.3828e- 

-01 

2.9602 

0.06 

0.6595e-03 

0.9864e- 

-01 

3.  1076 

0.08 

0.3260e-02 

0.1617 

3.2305 

0.  10 

0.7146e-02 

0.2272 

3.3368 

0.  12 

0. 1236e-01 

0.2946 

3.4303 

0.  14 

0. 1894e-01 

0.3638 

3.5125 

0.  16 

0.2693e-01 

0.4346 

3.5839 

0.18 

0.3634e-01 

0.5067 

3.6444 
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Wnrk-ftnergy  halanrp  r.hftrk: 

t  d(KE)  d{W) 

0.00  O.OOOOe+00  O.OOOOe+00 

0.02  0.8104e-02  0.9018e-02 

0.04  0.2120e-01  0.2326e-01 

0.06  0.3968e-01  0.4301e-01 

0.08  0.6401e-01  0.6873e-01 

0.10  0.9464e-01  0.1009 

0.12  0.1320  0.1398 

0.14  0.1765  0.1860 

0.16  0.2285  0.2397 

0.18  0.2884  0.3015 
where 

d(KE)  is  increase  in  kenetic  energy, 

and  d(W)  is  the  work  done  on  the  system. 
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(1-2)  Tnverse  dynamirs  problem: 

T n  1 1- i a  1 cnndl  tions  = 

q2=0.0000e+00:   q2=0.2094:   qj=0 . 5179e-01 ; 
q2=0.0000e  +  00:   q2=-0.1396;   q2="0-^^'^^: 
q3=0.0000e+00:   q3=-0 . 6981e-01 ;   q3=2 . 2970 : 

Generalized  coordinates,  velocities,  and  accelerations 
at  each  instant  are  specified  by  the  output  from  the 
forward  dynamic  problem. 


Resu 

Its: 

t 

Qi 

Q2 

Q3 

0. 

.00 

0. 

.6000e+01 

0. 

. 5000e+01 

0, 

. 4000e+01 

0. 

.02 

0, 

.7414e+01 

0, 

.6131e+01 

0, 

.4849e+01 

0, 

.04 

0, 

.8000e+01 

0, 

.  6600e  +  01 

0. 

.5200e+01 

0. 

.06 

0. 

.8449e+01 

0, 

.6960e+01 

0, 

.  5470e  +  01 

0. 

,08 

0. 

.8828e+01 

0. 

.7263e+01 

0. 

. 5697e+01 

0. 

10 

0, 

,9162e+01 

0. 

,7530e+01 

0. 

,5897e+01 

0. 

12 

0. 

,9464e+01 

0, 

.7771e+01 

0. 

,6078e+01 

0. 

14 

0. 

.9742e+01 

0. 

.7993e+01 

0. 

.6245e+01 

0. 

16 

0. 

. lOOOe+02 

0. 

.  8200e+01 

0, 

, 6400e+01 

0. 

18 

0. 

1024e+02 

0. 

8394e+01 

0. 

, 6546e+01 
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t 

Qi 

Q2 

Q3 

0. 

00 

0. 

6000E+01 

0. 

5000E+01 

0. 

4000e+01 

0. 

02 

0. 

7414e+01 

0. 

6131e+01 

0. 

4849e+01 

0. 

04 

0. 

8000e+01 

0. 

6600e+01 

0. 

.5200e+01 

0. 

06 

0, 

.8450e+01 

0. 

.6960e+01 

0. 

,  5470e+01 

0. 

,08 

0. 

. 8828e+01 

0. 

.  7263e+01 

0. 

. 5697e+01 

0. 

.  10 

0, 

.9162e+01 

0, 

.  7530e+01 

0. 

.5897e+01 

0, 

.  12 

0 

.9464e+01 

0 

.7771e+01 

0. 

.6078e+01 

0, 

.  14 

0 

. 9742e+01 

0 

. 7993e+01 

0 

.6245e+01 

0, 

.  16 

0 

.  lOOOe  +  02 

0 

. 8200e+01 

0 

.6400e+01 

0 

.18 

0 

. 1024e+02 

0 

.8394e+01 

0 

.6546e+01 

45 


rigora  1 
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F.XAMPI.F.  9:  4-T.TNK  ROBOT 

The  initial  position  or  home  position  and  coordinate 
systems  are  shown  in  figrue  2.  The  only  load  is  the  self- 
weight  of  the  members;  the  simulation  time  is  from  0.00  to 
0.18  second . 

The  input  data  is  given  below,  and  all  units  are  in 
MKS. 

Lengths : 

Li=1.5;   L2=1.0;   L3=0.75:   L4=0.25; 
Masses : 

mj=50;   m2=30;   m3=20;   m4=15; 
Mnments  nf  ineT-tia: 

Ilj^=9.375:    Iiy=1.0;    1-^^=9.31^; 

l2x=0.15:    l2y=2.5:    Ii^=2.5: 

l3x=01:    l3y=0.9375:    13^=0.9375; 

l4x=0.075;    l4y=0.8;    l4z=0.8; 
Centers  nf  mass: 

ri^=0;   riy=0.75;   ri^=0; 

r2x=0.50:   r2y=0;   v^^^O; 

r3^=0.375:   r3y=0;   r3^=0; 

^4x=-015:   r4y=0;   r42:=0: 
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(2-1)  Forward  Hynamir  problem: 

Tni  tial cnndi  tinns : 

qj=0;        q2=0:        ^3=0:       ^4=0: 

Ql=0.3491:       q2=0-2618;       q3=-0.2618;       44=0.08727 

GenerallT-ftd     fnrr.p    fimrt.ion    at    each — instant  : 

Qj=     10.0»*t^^2+6.0; 

Q2=    8.0*«t^/2+5  0; 

Q3=    6.0»*t^^2+4.0: 

Q^=    4.0»tt^''2+3.0; 
where    t    is    the    simulation    time. 


Rfisu 

Us: 
For  linV  1 : 

t 

^1 

^1 

0.00 

O.OOOOe+00 

0.3491 

0.02 

-0.6928e-02 

0.3436 

0.04 

0. 1374e-01 

0.3376 

0.06 

0.2043e-01 

0.3314 

0.08 

0.2699e-01 

0.3249 

0.  10 

0.3343e-01 

0.3182 

0.  12 

0.3972e-01 

0.3114 

0.14 

0.4588e-01 

0.3045 

0.  16 

0.5190e-01 

0.2974 

0.  18 

0.5778e-01 

0.2903 

^1 

-0.2486 
-0.2920 
-0.3086 
-0.3206 
-0.3302 
-0.3382 
-0.3450 
-0.3508 
-0.3560 
-0.3604 
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Fnr  link  2: 

t 

^2 

^2 

0.00 

O.OOOOe+00 

0. 

,2618 

0.02 

0.5246e-02 

0. 

,2628 

0.04 

0. 1051e-01 

0. 

2640 

0.06 

0. 1580e-01 

0. 

,2652 

0.08 

0.2112e-01 

0. 

.2664 

0.  10 

0.2646e-01 

0, 

,2676 

0.  12 

0.3183e-01 

0 

.2689 

0.  14 

0.3722e-01 

0 

.2702 

0.  16 

0.4263e-01 

0 

.2715 

0.  18 

0.4808e-01 

0 

.2727 

"^2 
0.4438e-01 

0.5564e-01 

0.5920e-01 

0.6136e-01 

0.6275e-01 

0.6361e-01 

0.6405e-01 

0.6417e-01 

0.6400e-01 

0.6361e-01 


Fnr  link  3: 

t 

"^3 

^3 

•  • 

^3 

0.00 

O.OOOOe+00 

-0.2618 

.0492 

0.02 

-0.5009e-02 

-0.2386 

.2471 

0.04 

-0.9527e-02 

-0.2130 

.3221 

0.06 

-0. 1352e-01 

-0. 1862 

.3749 

0.08 

-0.1697e-01 

-0.1584 

.4151 

0.  10 

-0. 1985e-01 

-0. 1299 

.4466 

0.  12 

-0.2216e-01 

-0. 1008 

.4711 

0.  14 

-0.2388e-01 

-0.7130e- 

-01 

.4900 

0.  16 

-0.2501e-01 

-0.4144e- 

-01 

,5038 

0.  18 

-0.2554e-01 

-0. 1135 

,5130 
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Fnr  link  4= 

t 

^4 

^4 

^'4 

0.00 

O.OOOOe+00 

0.8727e- 

-01 

0.4361 

0.02 

0. 1836e-02 

0.9639e- 

-01 

0.4716 

0.04 

0.3858e-02 

0.1059 

0.4854 

0.06 

0.6075e-02 

0.1157 

0.4964 

0.08 

0.8489e-02 

0. 1257 

0.5062 

0.  10 

O.lllOe-01 

0.1359 

0.5157 

0.  12 

0.1393e-01 

0. 1463 

0.5253 

0.  14 

0.1696e-01 

0.1569 

0.5352 

0.  16 

0.2020e-01 

0.1677 

0.5457 

0.  18 

0.2367e-01 

0. 1787 

0.5568 

Wnrk-energy  balance  check: 

t  d(KE) 

0.00  O.OOOOe+00 

0.02  0.6028e-01 

0.04  0.1306 

0.06  0.2083 

0.08  0.2930 

0.10  0.3842 

0.12  0.4819 

0.14  0.5861 

0.16  0.6967 

0.18  0.8137 


d(W) 
O.OOOOe+00 
0.6137e-01 
0.1329 
0.2119 
0.2977 
0.3901 
0.4890 
0.5944 
0.7062 
0.8244 
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where 

d(KE)  is  increase  in  kenetic  energy, 
and  d(W)  is  the  work  done  on  the  system. 
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(2-2)  TnvRrsft  dynamics  problem: 
Initial rnndi  t i  nns : 


qi=0.3491:    qj=-0.2486; 

q2=0.2618;    q2=0. 04438; 

q3=:-0.2618;    q3=1.0492; 

q4=0. 08727;    q4=0.4361 
Generalized  coordinates,  velocities,  and  accelerations 
at  each  instant  are  specified  by  the  output  from  the 
forward  dynamic  problem. 


qi=0.0000 
q2=0.0000 
q3=0.0000 
q4=0.0000 


Resu 1 t  s : 
t         Qi 

0.00  0.6000e+01 

0.02  0.7414e+01 

0.04  0.7999e+01 

0.06  0.8448e+01 

0.08  0.8827e+01 

0.10  0.9160e+01 

0.12  0.9462e+01 

0.14  0.9739e+01 

0.16  0.9997e+01 

0.18  0.1024e+02 


Q2 
0.5000e+01 

0.6131e+01 

0.6600e+01 

0.6960e+01 

0.7263e+01 

0.7530e+01 

0.7772e+01 

0.7994e+01 

0.8201e+01 

0.8395e+01 


Q3 

0.4000e+01 

0.4848e+01 
0.5200e+01 
0.5469e+01 
0.5696e+01 
0.5897e+01 
0.6078e+01 
0.6244e+01 
0.6399e+01 
0.6545e+01 


Q4 

0.3000e+01 

0.3566e+01 
0.3800e+01 
0.3980e+01 
0.4131e+01 
0.4265e+01 
0.4386e+01 
0.4497e+01 
0.4600e+01 
0.4697e+01 
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Thf.    ftvart  values  as  givpn  by T.hfi  functions 


0.00  0.6000E+01 

0.02  0.7414e+01 

0.04  0.8000e+01 

0.06  0.8450e+01 

0.08  0.8828e+01 

0.10  0.9162e+01 

0. 12  0.9464e+01 

0.14  0.9742e+01 

0.16  O.lOOOe+02 

0.18  0.1024e+02 


Q2 

0.5000E+01 

0.6131e+01 
0.6600e+01 
0.6960e+01 
0.7263e+01 
0.7530e+01 
0.7771e+01 
0.7993e+01 
0.8200e+01 
0.8394e+01 


Q3 

0.4000e+01 

0.4849e+01 
0.5200e+01 
0.5470e+01 
0.5697e+01 
0.5897e+01 
0.6078e+01 
0.6245e+01 
0.6400e+01 
0.6546e+01 


Q4 

0.3000e+01 

0.3566e+01 
0.3800e+01 
0.3980e+01 
0.4131e+01 
0.4265e+01 
0.4386e+01 
0.4497e+01 
0.4600e+01 
0.4697e+01 
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CHAPTER    6 
CONCLUSIONS 

The  major  objective  of  this  research  was  to  develop  an 
automated  simulator  generator  for  general  serial-link  open- 
chain  robots  by  using  the  symbolic  language  REDUCE-3.  The 
simulator,  based  on  the  Lagrangian  formulation,  can  solve 
both  the  forward  and  the  inverse  dynamic  problems  of  any 
robot  specified  by  the  user,  provided  only  revolute  and 
translat ional     joints    are    used. 

The  system  equations  are  set  up  for  user-specified 
robot  and  solved  without  any  modifications  or 
simplifications  so  that  exact  simulations  (within  round- 
off) can  be  obtained,  the  solution  of  the  system  equations 
is  done  in  a  way  that  avoids  excessively  large  explicit 
expressions  and  greatly  reduces  the  computing  time.  A  more 
efficient  way  of  using  symbolic  languages  was  also 
developed  through  the  compromise  between  REDUCE  and  FORTRAN 
languages,     as    suggested    in    [18]. 

Several  advantages  come  with  the  use  of  symbolic 
processing  in  this  application.  First  of  all,  it 
eliminates  the  time-consuming  and  tedious  manual  derivation 
process.  Secondly,  the  simulators  generated  would  be  very 
reliable     since     the     possibility     of     programmer     error     is 
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minimized.  Another  benefit  is  the  remarkable  flexibility 
of  the  simulator  generation  program,  because  this  program 
computes  all  quantities  symbolically  rather  than 
numerically,  it  can  handle  a  wide  variety  of  problems 
without  modification.  This  become  important  when  attempting 
to  extend  the  scope  of  the  simulators  to  include  controls 
and    optimization. 

The  work  presented  in  this  thesis  offers  wide  scope  for 
future  development.  The  techniques  used  to  express  and 
solve  the  Lagrangian  equation  can  be  applied  to  various 
mechanical  systems.  since  most  real-world  mechanical 
dynamics    problem    can    be    described    by    the    Lagrangian. 

The  method  for  describing  the  geometry  of  robots  can  be 
improved  and  standardized  by  adopting  the  widely  used 
Denav i t-Har t enbe rg      convention.  The      schemes      used      for 

calculating  complicated  terms  in  the  system  equation  could 
be  further  improved  by  treating  the  matrices  in  a  more 
efficient  way  and  by  finding  a  better  trade-off  between 
computing  time  and  memory  space.  Due  to  the  presumption 
that  only  one  degree  of  freedom,  or  one  generalized 
coordinate  is  associated  with  each  link,  the  simulators 
generated  becomes  clumsy  when  simulating  a  multi-degree  of 
freedom  link,  in  which  case  a  single  link  must  be  treated 
as  several  artificial  links.  Improvements  on  this  defect 
can  result  in  better,  more  efficient  simulation  of  wrist 
joints    and    other    higher    order    pairs. 
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In  the  context  of  on-line  dynamic  control,  it  was 
suggested  in  [19]  that  the  canonical  formulation  has 
advantages  over  Lagrange's  equations  of  motion,  since  it 
avoids  taking  time  derivatives,  which  in  turn  results  in 
simpler  expressions  in  the  equations  of  motion.  Also,  the 
use  of  super  computers,  array  processors,  or  parallel 
processing  computers  would  help  in  obtaining  numeric 
solutions  within  the  sampling  time. 

There  is  ample  scope  for  development  in  computer- 
related  areas  as  well.  Computer  graphics  can  be  used  to 
display  the  results  of  dynamic  analysis.  Interactive 
commands  could  also  be  developed  to  make  the  program  more 
user- f r i end  1 y  and  to  give  the  user  more  control  over  the 
design  process.  An  expert  system  for  design  purposes  could 
be  included  in  the  preprocessor  to  determine  the  basic 
profile  and  mechanical  requirements  of  robots  according  to 
their  uses. 

Since  it  is  nearly  impossible  to  standardize  problems 
in  the  design  of  mechanical  systems.  Automatic  simulator 
generators  are  extremely  attractive  as  tools  for  the 
simulation,  analysis,  design  ,  and  optimization  of  robot 
manipulators;  they  are  also  a  valuable  aid  in  planning  the 
efficient  utilization  of  these  manipulators  for  specific 
tasks.  This  is  certainly  an  area  worthy  of  more  research 
work  and  the  one  which  will  become  increasingly  popular  in 
the  near  future. 
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APPENDIX  A 
The  input  to  the  REDUCE  based  simulator  generator  for 
both  examples  are  given  below. 

For  the  3-R  ROBOT  in  example  1: 

NUML:=3: 

L(1):=RL(1);  L(2):=RL(2):  L(3):=RL(3); 
CM{1 , 1):=XCM(1) :  CM(1.2):=0:  CM(1,3):=0:  CM(1.4)=1: 
CM{2. 1) :=XCM(2);  CM(2,2):=0;  CM(2.3):=0:  CM(2,4)=1: 
CM(3. 1) :=XCM(3);  CM(3.2):=0:  CM(3,3):=0:  CM(3.4)=1; 
MAS{1):=RM{1) :  MAS(2) : =RM(2) ;  MAS{3) : =RM(3) ; 
MI{1.1):=RMI(1. 1):  MI{1.2):=RMI(1.2):  MI ( 1 . 3) : =RMI ( 1 . 3) 
MI(2, 1):=RMI(2. 1) :  MI (2 . 2) : =RMI (2 . 2) :  MI (2 , 3) : =RMI (2 . 3) 
MI(3. 1):=RMI{3, 1) ;  MI (3 . 2) : =RMI (3 . 2) :  MI (3 . 3) : =RMI (3 , 3 ) 
DIRECT(l) :=X:  DIRECT( 2) : =X :  DIRECT{3) : =X: 
JOINT(l) :=REV:  J0INT(2) : =REV ;  J0INT(3) : =REV ; 
AXIS{1):=Z;  AXIS{2):=Z;  AXIS{3):=Z:  -.END; 
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For  the  4-LINK  ROBOT  In  example  2: 


NUML:=4: 

L{1):=RL(1):  L(2):=RL(2);  L(3):=RL(3):  L(4):=RL(4): 


CM(1.1):=0;  CM( 1 . 2) : =YCM{ 1 ) ;  CM(1.3 
CM(2.1):=XCM{2) :  CM(2.2):=0;  CM(2.3 
CM(3.1):=XCM(3) :  CM(3.2):=0:  CM(3.3 
CM(4.1):=XCM(4):  CM(4.2):=0;  CM{4.3 
MAS(l) :=RM{1):  MAS( 2) : =RM(2) ;  MAS{3 


MI{1 . 1):=RMI{1. 1) 
MI(2. 1):=RMI(2. 1) 
MI(3. 1) :=RMI(3. 1) 
MI(4. 1) :=RMI(4,1) 


MI(1.2) :=RMI(1.2 
MI(2.2):=RMI(2.2 
MI(3.2) :=RMI(3,2 


:=0:  CM{1.4)=1: 
:=0:  CM(2.4)=1: 
:=0:  CM(3.4)=1: 
:=0:  CM(4.4)=1; 
:=RM(3);  MAS(4) : =RM(4) ; 
;  MI(1.3) :=RMI{1 ,3) 
:  MI(2.3) :=RMI(2.3) 
;  MI(3.3):=RMI(3.3) 
;  MI{4.3) :=RMI(4.3) 


MI(4.2):=RMI(4.2; 
DIRECT(1):=Y;  DIRECT(2) : =X ;  DIRECT{3) : =X ;  DIRECT(4) : =X ; 
J0INT{1) :=REV;  J0INT(2) : =REV :  J0INT(3) : =REV ;  J0INT(4) : =TRAN 
AXIS(1):=Y:  AXIS(2):=Z:  AXIS{3):=Y;  ;END; 
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APPENDIX  B 
The  meaning  of  all  the  variables  in  APPENDIX  A  is  given 
in  this  section  as  follows. 
NUML  :  the  number  of  links 
L     :  the  length  of  the  link 

CM    :  the  location  of  the  center  of  mass  of  the  link 
MAS   :  the  mass  of  the  link 
MI    :  the  moment  of  inertia  of  the  link 
DIRECT  :  the  direction  in  which  the  link  is  pointing 
JOINT   :  the  joint  type 
AXIS    :  the  axis  of  rotation  of  the  link 
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ABSTRACT 

This  thesis  presents  a  generalized  technique  for 
generating  simulators  for  serial,  open-chain  robot 
manipulators  using  symbolic  computing.  These  simulators 
can  handle  both  direct  dynamics  problem  as  well  as  the 
inverse  dynamics  problem.  In  both  cases.  the  general 
governing  equations  are  derived  using  joint  coordinates  and 
the  principles  of  Lagrangian  mechanics.  These  equations 
are  then  recast  into  a  form  that  is  better  suited  for 
implementation  in  the  symbolic  processing  language  REDUCE. 
A  symbolic  processor,  written  in  REDUCE,  has  been  set  up  to 
evaluate  all  the  required  expressions  in  these  equations 
and  to  write  these  out  in  the  form  of  a  complete  FORTRAN 
simulation  program;  this  FORTRAN  simulator  is  then  executed 
to  obtain  the  system  response.  Numerical  examples  are 
presented  to  demonstrate  the  efficiency  of  the  method  and 
to  establish  the  feasibility  of  the  proposed  technique  of 
generating  simulators  automatically. 

The  detail  of  design  sensitivity  analysis  is  developed 
so  that  it  is  ready  for  implementing  to  the  symbolic 
computing  language.  Combining  the  sensitivity  and  dynamics 
analysis,  the  simulator  generated  is  then  able  to  be 
extended  to  include  the  design  optimization/update  part  for 
improved  designs. 


